{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# B 08-7 随机区组设计的方差分析 Two way ANOWA\n",
    "\n",
    "## 案例\n",
    "\n",
    "将大鼠分成 10 各区组，同一区组的大鼠随机分进 3 个实验组，分别予以不同处理，一定时间后检测睾丸 MT 含量。比较三种不同处理的大鼠 MT 有无区别？\n",
    "\n",
    "### 数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id</th>\n",
       "      <th>house</th>\n",
       "      <th>group</th>\n",
       "      <th>value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>control</td>\n",
       "      <td>40.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>control</td>\n",
       "      <td>44.8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>control</td>\n",
       "      <td>36.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>control</td>\n",
       "      <td>49.9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>control</td>\n",
       "      <td>59.8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   id  house    group  value\n",
       "0   1      1  control   40.6\n",
       "1   2      2  control   44.8\n",
       "2   3      3  control   36.7\n",
       "3   4      4  control   49.9\n",
       "4   5      5  control   59.8"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.formula.api import ols\n",
    "\n",
    "df = pd.read_csv(\"B_08_7.csv\")\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLK0lEQVR4nO3deXyMd+IH8M9kcorMREKOIbcggjiCBm1V0tKWpbWuTQlRtl2WuGv7Cy11VlB11VFHUbVtqdotTUNdjYhINEgJgiBHK8lMDjlnfn9Y006TkInMPDNPPu/Xa15rvs8zz3wetnk+eU6JRqPRgIiIiEikLIQOQERERGRILDtEREQkaiw7REREJGosO0RERCRqLDtEREQkaiw7REREJGosO0RERCRqlkIHMAVqtRr37t2Dg4MDJBKJ0HGIiIioDjQaDQoLC6FQKGBhUfv+G5YdAPfu3YOHh4fQMYiIiKgeMjMz0apVq1qns+wAcHBwAPDwL0smkwmchoiIiOpCpVLBw8NDux2vDcsOoD10JZPJWHaIiIjMzJNOQeEJykRERCRqLDtEREQkaiw7REREJGo8Z6eO1Go1ysvLhY7RaFlZWUEqlQodg4iIzBDLTh2Ul5cjIyMDarVa6CiNmqOjI9zc3HgvJCIi0gvLzhNoNBpkZWVBKpXCw8PjsTctIsPQaDQoKSlBbm4uAMDd3V3gREREZE5Ydp6gsrISJSUlUCgUaNKkidBxGi07OzsAQG5uLlxcXHhIi4iI6oy7KZ6gqqoKAGBtbS1wEnpUNisqKgROQkRE5oRlp454nojw+G9ARET1wbJDREREosZzdqhRuF5wHd/f/B4aaPCi14vwb+YvdCQiIjISlh0SvZ2XdmLFuRXQQAMA2HBhA6Z2nYo3O74pcDIiIjIGHsYiUcsqysLKpJXaovPImvNrcFt1W6BURERkTCw7IsO7POv68c6PqNJUVRvXQINjmccESERERMbGsmPiCgsLER4eDnt7e7i7u2PVqlXo27cvoqKiAADe3t5YuHAhxowZA5lMhokTJwIAvvrqKwQGBsLGxgbe3t6IiYnRWa5EIsGBAwd0xhwdHbF9+3YAwM2bNyGRSLB371706tULtra26NChA44fP27oVW5Q1ha13zLAWsrbCRARNQYsOyZu+vTpOH36NA4ePIjY2FicPHkS58+f15lnxYoVCAoKQnJyMqKjo5GUlIThw4dj5MiRSE1NxXvvvYfo6GhtkdHHrFmzMGPGDCQnJyMkJASDBg3C/fv3G2jtDC/UMxS2Uttq41YWVnjR60UBEhERkbGx7JiwwsJC7NixAytWrEBoaCg6dOiAbdu2aW90+Ei/fv0wY8YM+Pn5wc/PDytXrkRoaCiio6PRpk0bjB07FpMnT8aHH36od4bJkydj6NChCAgIwIYNGyCXy7F169aGWkWDc7R1xNLnlqKJ5e93v7aztMOSZ5eguV1zAZMREZGx8GosE3bjxg1UVFSgR48e2jG5XI62bdvqzBccHKzzPi0tDYMHD9YZ6927N1avXo2qqiq9HrUQEhKi/bOlpSWCg4ORlpamz2oILtQzFEeHH8XJuyeh0WjQp2UfOFg7CB2LiIiMhGVHBOzt7fX+jEQigUaje4WSmB/DYG9ljwHeA4SOQUREAuBhLBPm6+sLKysrJCYmaseUSiWuXr362M8FBATg9OnTOmOnT59GmzZttHt1WrRogaysLO309PR0lJSUVFvWmTNntH+urKxEUlISAgIC6rU+REREQuCeHRPm4OCAiIgIzJo1C05OTnBxccH8+fNhYWHx2OdEzZgxA927d8fChQsxYsQIxMfHY+3atVi/fr12nn79+mHt2rUICQlBVVUV5syZAysrq2rLWrduHfz9/REQEIBVq1YhPz8fkZGRBllfIiIiQ+CeHRO3cuVKhISEYODAgQgLC0Pv3r0REBAAW9vqVxg90rVrV+zbtw979+5Fhw4dMG/ePCxYsABjx47VzhMTEwMPDw88++yz+Nvf/oaZM2dqnyr+R0uXLsXSpUsRFBSEU6dO4eDBg2jenCf2EhGR+ZBo/nziRiOkUqkgl8uhVCohk8l0ppWWliIjIwM+Pj6PLRjGUlxcjJYtWyImJgbjx4832PfcvHkTPj4+SE5ORufOnQ32PfowtX8LIiIS1uO233/Ew1gmLjk5Gb/88gt69OgBpVKJBQsWAEC1q62IiIioZiw7ZmDFihW4cuUKrK2t0a1bN5w8eZKHkoiIiOqIZcfEdenSBUlJSUb/Xm9v72qXphMREZkjnqBMREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskMmYfv27XB0dBQ6BhERiRDLDtWbt7c3Vq9eLXQMIiKix2LZIYOqqqqCWq0WOgYRETViLDtGcuRSNgavPYWA6MMYvPYUjlzKNvh3qtVqLF++HK1bt4aNjQ08PT2xaNEiAEBqair69esHOzs7ODs7Y+LEiSgqKtJ+duzYsRgyZAhWrFgBd3d3ODs7Y9KkSaioqAAA9O3bF7du3cK0adMgkUi0T2F/dDjq4MGDaN++PWxsbHD79m3k5+djzJgxaNasGZo0aYKXX34Z6enpBv87ICIiYtkxgiOXsvH3z5Jw4Y4SDyqqcOGOEm/tSjJ44Zk7dy6WLl2K6OhoXL58GXv27IGrqyuKi4vRv39/NGvWDImJifj3v/+NH374AZMnT9b5/LFjx3D9+nUcO3YMO3bswPbt27F9+3YAwNdff41WrVphwYIFyMrKQlZWlvZzJSUlWLZsGbZs2YJLly7BxcUFY8eOxblz53Dw4EHEx8dDo9HglVde0ZYnIiIiQ+HjIoxg/bFr1cY0GmD9j9fRP9DNIN9ZWFiIjz76CGvXrkVERAQAwM/PD3369MHmzZtRWlqKnTt3wt7eHgCwdu1aDBo0CMuWLYOrqysAoFmzZli7di2kUinatWuHV199FXFxcZgwYQKcnJwglUrh4OAANzfddaioqMD69esRFBQEAEhPT8fBgwdx+vRp9OrVCwCwe/dueHh44MCBAxg2bJhB/g6IiIgAgffsnDhxAoMGDYJCoYBEIsGBAwdqnfett96CRCKpdkJsXl4ewsPDIZPJ4OjoiPHjx+scjjEFV3NqzpOeU2iw70xLS0NZWRlCQ0NrnBYUFKQtOgDQu3dvqNVqXLlyRTsWGBgIqVSqfe/u7o7c3Nwnfre1tTU6deqk832Wlpbo2bOndszZ2Rlt27ZFWlqa3utGRESkD0HLTnFxMYKCgrBu3brHzrd//36cOXMGCoWi2rTw8HBcunQJsbGxOHToEE6cOIGJEycaKnK9tHFtWuO4v6uDwb7Tzs7uqZdhZWWl814ikdTpZGM7OzvtOTxERERCE7TsvPzyy/jggw/w2muv1TrP3bt38c9//hO7d++utvFNS0vD4cOHsWXLFvTs2RN9+vTBxx9/jL179+LevXuGjl9n/3ihNf687ZdIgEl9/Qz2nf7+/rCzs0NcXFy1aQEBAbhw4QKKi4u1Y6dPn4aFhQXatm1b5++wtrZGVVXVE+cLCAhAZWUlEhIStGP379/HlStX0L59+zp/HxERUX2Y9AnKarUao0ePxqxZsxAYGFhtenx8PBwdHREcHKwdCwsLg4WFhc6G9c/KysqgUql0XobUP9ANG9/ohiAPRzSxliLIwxGfvNENLxnofB0AsLW1xZw5czB79mzs3LkT169fx5kzZ7B161aEh4fD1tYWERERuHjxIo4dO4Z//vOfGD16tPZ8nbrw9vbGiRMncPfuXfz222+1zufv74/BgwdjwoQJOHXqFC5cuIA33ngDLVu2xODBgxtidYmIiGpl0icoL1u2DJaWlpgyZUqN07Ozs+Hi4qIzZmlpCScnJ2Rn136l05IlS/D+++83aNYn6R/oZrCTkWsTHR0NS0tLzJs3D/fu3YO7uzveeustNGnSBEeOHMHUqVPRvXt3NGnSBEOHDsXKlSv1Wv6CBQvw97//HX5+figrK4NGo6l13m3btmHq1KkYOHAgysvL8dxzz+G///1vtb11REREDU2iedwWyogkEgn279+PIUOGAACSkpLw6quv4vz589pzdby9vREVFYWoqCgAwOLFi7Fjxw6dk2oBwMXFBe+//z7efvvtGr+rrKwMZWVl2vcqlQoeHh5QKpWQyWQ685aWliIjIwM+Pj6wtbVtoLWl+uC/BRER/ZFKpYJcLq9x+/1HJnsY6+TJk8jNzYWnpycsLS1haWmJW7duYcaMGfD29gYAuLm5Vbs6qLKyEnl5edUuh/4jGxsbyGQynRcRERGJk8kexho9ejTCwsJ0xvr374/Ro0dj3LhxAICQkBAUFBQgKSkJ3bp1AwAcPXoUarVa5zJnIiIiarwELTtFRUW4du33G+5lZGQgJSUFTk5O8PT0hLOzs878VlZWcHNz014xFBAQgAEDBmDChAnYuHEjKioqMHnyZIwcObLGy9SJiIio8RH0MNa5c+fQpUsXdOnSBQAwffp0dOnSBfPmzavzMnbv3o127dohNDQUr7zyCvr06YNNmzYZKjIRERGZGUH37PTt2/exV/D82c2bN6uNOTk5Yc+ePQ2YioiIiMTEZE9QJiIiImoILDtEREQkaiw7REREJGosO0RERCRqLDtEREQkaiw7JJgpU6agW7dusLGxQefOnYWOQ0REIsWyQ1rl5eVG/87IyEiMGDHC6N9LRESNB8uOsaQdAja9ACxyf/i/aYcM/pWFhYUIDw+Hvb093N3dsWrVKvTt21f7IFVvb28sXLgQY8aMgUwmw8SJEwEAX331FQIDA2FjYwNvb2/ExMToLFcikeDAgQM6Y46Ojti+fTuAh/dDkkgk2Lt3L3r16gVbW1t06NABx48f1/nMmjVrMGnSJPj6+hpk/cm0xd2Ow6hDo9Bjdw/87T9/w9HbR4WOREQixbJjDGmHgC/CgXvngYqSh//7xRsGLzzTp0/H6dOncfDgQcTGxuLkyZM4f/68zjwrVqxAUFAQkpOTER0djaSkJAwfPhwjR45Eamoq3nvvPURHR2uLjD5mzZqFGTNmIDk5GSEhIRg0aBDu37/fQGtH5izudhyijkXh4v2LeFD5AKm/pSLqWBSO3T4mdDQiEiGWHWM4GVPDoAY4tdJgX1lYWIgdO3ZgxYoVCA0NRYcOHbBt2zZUVVXpzNevXz/MmDEDfn5+8PPzw8qVKxEaGoro6Gi0adMGY8eOxeTJk/Hhhx/qnWHy5MkYOnQoAgICsGHDBsjlcmzdurWhVpHM2OafN1cb00CDLalbBEhDRGLHsmMMv/5S83huLeMN4MaNG6ioqECPHj20Y3K5XPsQ1UeCg4N13qelpaF37946Y71790Z6enq1ovQkISEh2j9bWloiODgYaWlpei2DxOl6wfUax68VXKtx3NQUlBZgS+oWzDw+E2vOr0FWUZbQkYjoMQR9Nlaj0aLdw0NXf+bSzvhZ/sTe3l7vz0gkkmrPNKuoqGioSNQI+Dn64dL9S9XGWzu2FiCNfrKKsjD6u9HIKcnRjn3+y+fY0n8LAp0DBUxGRLXhnh1jeHYGAMmfBiX/GzcMX19fWFlZITExUTumVCpx9erVx34uICAAp0+f1hk7ffo02rRpA6lUCgBo0aIFsrJ+/002PT0dJSUl1ZZ15swZ7Z8rKyuRlJSEgICAeq0PicuEjhOqjUkgwfiO4wVIo59Pfv5Ep+gAQFFFEVaeM9xhaSJ6OtyzYwwBA4ERux6eo5P7y8M9Os/OANq9arCvdHBwQEREBGbNmgUnJye4uLhg/vz5sLCwgETy5+L1uxkzZqB79+5YuHAhRowYgfj4eKxduxbr16/XztOvXz+sXbsWISEhqKqqwpw5c2BlZVVtWevWrYO/vz8CAgKwatUq5OfnIzIyUjv92rVrKCoqQnZ2Nh48eICUlBQAQPv27WFtbd1wfxlkckK9QrG672psTt2MG8ob8JP7YUKnCejn2U/oaE8Ufy++xvGz2WdRoa6AlUX1/xaISFgsO8YSMPDhy4hWrlyJt956CwMHDoRMJsPs2bORmZkJW1vbWj/TtWtX7Nu3D/PmzcPChQvh7u6OBQsWYOzYsdp5YmJiMG7cODz77LNQKBT46KOPkJSUVG1ZS5cuxdKlS5GSkoLWrVvj4MGDaN68uXb6m2++qXM5epcuXQAAGRkZ8Pb2fvq/ADJpoV6hCPUKFTqG3mQ2Mtwrvldt3MHKAZYS/kglMkUSzZ9PvmiEVCoV5HI5lEolZDKZzrTS0lJkZGTAx8fnsSXBHBQXF6Nly5aIiYnB+PGGO1xw8+ZN+Pj4IDk5uUHvjCymfwsyX5//8jkWJyyuNj6m/RjM6j5LgEREjdfjtt9/xF9DRCw5ORm//PILevToAaVSiQULFgAABg8eLHAyIvM1su1I3Cm8g89/+RwV6gpYSCzQ37s/pnSdInQ0IqoFy47IrVixAleuXIG1tTW6deuGkydP6hxKIiL9SCQSzOo+C+M7jsf1guvwcPCAm72b0LGI6DFYdkSsS5cuNZ5LY2je3t7VLk0nEhsnWyc4uTkJHYOI6oCXnhMREZGosewQERGRqLHsEBERkaix7BCZCWWZEqpyldAxiIjMDk9QJjJxN5Q3sOjMIpzNPgsJJOjdsjfmPTMP7k3dhY5GRGQWuGeHyIQ9qHyACUcm4Gz2WQCABhqcunsKE2Mnokqt31Poicj8XMm7guk/TkfovlCMOjQK317/VuhIZol7dohM2Pc3v0fug9xq4zdVN3Hq7ik87/G8AKmIyBgylBmIOByB4opiAEDug1z869S/oCxT4o32bwiczrxwzw4J4sKFCxg1ahQ8PDxgZ2eHgIAAfPTRR0LHMjlZxVn1mkZE5m/HpR3aovNHm1M3o0JdIUAi88U9O6RVXl5utKeNJyUlwcXFBbt27YKHhwd++uknTJw4EVKpFJMnTzZKBnPQqXmnWqd1bNHRiEmIyNiu5l+tcTyvNA+/lvwKRVOFkROZL+7ZMZK423EYdWgUeuzugVGHRiHudpzBv7OwsBDh4eGwt7eHu7s7Vq1ahb59+yIqKgrAwzsdL1y4EGPGjIFMJsPEiRMBAF999RUCAwNhY2MDb29vxMTE6CxXIpHgwIEDOmOOjo7Yvn07gIcPApVIJNi7dy969eoFW1tbdOjQQecJ55GRkfjoo4/w/PPPw9fXF2+88QbGjRuHr7/+2mB/H+YoRBGCnu49q42/6PUiAp0DBUhERMbiKfOscdzBygHOds5GTmPeWHaMIO52HKKOReHi/Yt4UPkAF+9fxLRj0wxeeKZPn47Tp0/j4MGDiI2NxcmTJ3H+/HmdeVasWIGgoCAkJycjOjoaSUlJGD58OEaOHInU1FS89957iI6O1hYZfcyaNQszZsxAcnIyQkJCMGjQINy/f7/W+ZVKJZycePv9P5JIJFgXug7Tuk1Dx+Yd0blFZ7zT4x0sf2650NGIyMBGB4yGpUX1AzCjAkbBRmojQCLzxcNYRrDl5y3VxjTQYGvqVoR6hhrkOwsLC7Fjxw7s2bMHoaEPv2Pbtm1QKHR3e/br1w8zZszQvg8PD0doaCiio6MBAG3atMHly5fx4YcfYuzYsXplmDx5MoYOHQoA2LBhAw4fPoytW7di9uzZ1eb96aef8MUXX+A///mPXt/RGNhIbRDZIRKRHSKFjkJERhTYPBAbwjbg4/Mf4+fffkYLuxb4W8Df+LOgHlh2jOC68nqN49cKrhnsO2/cuIGKigr06NFDOyaXy9G2bVud+YKDg3Xep6WlYfDgwTpjvXv3xurVq1FVVQWpVFrnDCEhIdo/W1paIjg4GGlpadXmu3jxIgYPHoz58+fjpZdeqvPyiYjE7hn3Z/DMq8+gSl0FqUXdf/6SLh7GMgI/uV+N460dWxs5SXX29vZ6f0YikVR7qnlFRf2uDLh8+TJCQ0MxceJE/N///V+9lkFEJHYsOk+HZccI3uz0JiSQ6IxJIMGbHd802Hf6+vrCysoKiYmJ2jGlUomrV2s+u/+RgIAAnD59Wmfs9OnTaNOmjXavTosWLZCV9ftlz+np6SgpKam2rDNnzmj/XFlZiaSkJAQEBGjHLl26hBdeeAERERFYtGiRfitIRERURzyMZQShnqFY9cIqbE3dimsF19DasTXe7Pgm+nn2M9h3Ojg4ICIiArNmzYKTkxNcXFwwf/58WFhYQCKR1Pq5GTNmoHv37li4cCFGjBiB+Ph4rF27FuvXr9fO069fP6xduxYhISGoqqrCnDlzYGVlVW1Z69atg7+/PwICArBq1Srk5+cjMvLhseaLFy+iX79+6N+/P6ZPn47s7GwAgFQqRYsWLRr4b4OIiBozlh0jCfUMNdjJyLVZuXIl3nrrLQwcOBAymQyzZ89GZmYmbG1ta/1M165dsW/fPsybNw8LFy6Ei5sL/jXvX4iIiNDOExMTg3HjxuHZZ5+FQqHARx99hKSkpGrLWrp0KZYuXYqUlBS0bt0aBw8eRPPmzQEAX375JX799Vfs2rULu3bt0n7Gy8sLN2/ebLi/BCIiavQkmj+ffNEIqVQqyOVyKJVKyGQynWmlpaXIyMiAj4/PY0uCOSguLkbLli0RExOD8ePH1zpfWWUZbqpuolJdqR2zs7SDl8yrTseNb968CR8fHyQnJ6Nz584NER2AuP4tiIjE5P6D+9h/bT8ylBnwd/THa/6vQW4jN/j3Pm77/UfcsyNiycnJ+OWXX9CjRw8olUosWLAAAKpdbfVnWcVZOkUHePhAyt9Kf4NrE1eD5SUiIvNzo+AGxh0Zh7zSPO3YZ2mfYceAHWjl0ErAZL/jCcoi9+imgWFhYSguLsbJkye1h5JqUqWuqvFZLABQWF5oqJhERGSmVp1fpVN0ACC3JBfrUtYJlKg67tkRsS5dutR4Ls3jSCQSSCCBBtWPbv75irLaeHt7V7s0nYiIxOmnuz/VOH767ukax4XAPTukw0JiAQcbhxqnGeP4KxERmZem1k31GhcCy04dNaY9Fe727rC11D0BWGYjg7OtsA+ea0z/BkRE5uK11q/pNS4EHsZ6gkc30isvL4ednZ3AaYzD0sISvnJflFSWoLyqHHaWdtXKjxAe3biwpnv6EBGRMP7R+R+4V3wPhzMOQwMNpBIp/uL3F4zrME7oaFq89ByPv3RNo9Hg9u3bqKiogEKhgIUFd4YZm0ajQUlJCXJzc+Ho6Ah3d3ehIxER0Z/cKbyDW6pb8HP0g5u9m1G+k5eeNxCJRAJ3d3dkZGTg1q1bQsdp1BwdHeHmZpz/gIiISD+tHFqZzKXmf8ayUwfW1tbw9/dHeXm50FEaLSsrK72euE5ERPQIy04dWVhY8K69REREZognoBAREZGosewQERGRqAladk6cOIFBgwZBoVBAIpHgwIED2mkVFRWYM2cOOnbsCHt7eygUCowZMwb37t3TWUZeXh7Cw8Mhk8ng6OiI8ePHo6ioyMhrQkRERKZK0LJTXFyMoKAgrFtX/fkZJSUlOH/+PKKjo3H+/Hl8/fXXuHLlCv7yl7/ozBceHo5Lly4hNjYWhw4dwokTJzBx4kRjrQIRERGZOJO5z45EIsH+/fsxZMiQWudJTExEjx49cOvWLXh6eiItLQ3t27dHYmIigoODAQCHDx/GK6+8gjt37kChUNTpu+t6nT4RERGZjrpuv83qnB2lUgmJRAJHR0cAQHx8PBwdHbVFBwDCwsJgYWGBhISEWpdTVlYGlUql8yIiIiJxMpuyU1paijlz5mDUqFHa9padnQ0XFxed+SwtLeHk5ITs7Oxal7VkyRLI5XLty8PDw6DZiYiISDhmUXYqKiowfPhwaDQabNiw4amXN3fuXCiVSu0rMzOzAVISERGRKTL5mwo+Kjq3bt3C0aNHdY7Jubm5ITc3V2f+yspK5OXlPfaxAjY2NrCxsTFYZiIiIjIdJr1n51HRSU9Pxw8//ABnZ2ed6SEhISgoKEBSUpJ27OjRo1Cr1ejZs6ex4xIREZEJEnTPTlFREa5du6Z9n5GRgZSUFDg5OcHd3R1//etfcf78eRw6dAhVVVXa83CcnJxgbW2NgIAADBgwABMmTMDGjRtRUVGByZMnY+TIkXW+EouIiIjETdBLz3/88Ue88MIL1cYjIiLw3nvvwcfHp8bPHTt2DH379gXw8KaCkydPxrfffgsLCwsMHToUa9asQdOmTeucg5eeExGZpwu/XsCuy7twp/AO2ju3R0RgBDxlnkLHIiOp6/bbZO6zIySWHSIi83PizglMPToVlZpK7ZiDtQN2vbILvnJfAZORsYjyPjtikZidiHdPvYupR6di35V9KKsqEzoSEZHZ+ej8RzpFBwAKywux+efNAiUiU2XyV2OJza7Lu7AscZn2/dHMo/hvxn+x+cXNsJJaCZiMiMh8lFSU4Gr+1RqnpeSmGDcMmTyWHSMqLC/EmuQ11caTcpJw+OZhDPIbJEAqIiJdKbkpOHLzCDTQ4EWvF9HNtZvQkaqxtbSFo40jCsoKqk1zs6/91iPUOPEwlhGl/pqKB5UPapyWkFX74y2IiIzl4+SPMfq70diVtgu703Zj7OGxWJ20WuhY1VhILDCi7Ygap4UHhBs5DZk6lh0jcrR1rHVaM9tmxgtCRFSDW6pbNZ7vsvXiVtwouCFAosd7O+htRHaIhL2VPQCguV1zRD8TjTCvMIGTkanhYSwjau/cHoHOgbh0/5LOuKXEEq+1fk2gVERED526ewoa1HyB7ok7J+DraFpXOEktpJjWbRreDnob+aX5aNGkBSwtuFmj6rhnx8hWv7AaXV26at+3sGuB5c8vN7kfIkTU+DSxbFL7NKvapwnN1tIW7k3dWXSoVrzPDoS5z84t1S0UlRehrVNb/gdKRCZBVa7CS1++hOKKYp1xO0s7xP41FnIbuUDJiGrG++yYOC+ZFwKbB7LoEJHJkFnLsKrvKjjZOmnHnGydsKrvKhYdMmvc0hIRkVaIIgQ//PUHJGYnQgMNurt1h7XUWuhYRE+FZYeIiHRYSa3Qq2UvoWMQNRgexiIiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiUWPZISIiIlFj2SEiIiJRY9khIiIiURO07Jw4cQKDBg2CQqGARCLBgQMHdKZrNBrMmzcP7u7usLOzQ1hYGNLT03XmycvLQ3h4OGQyGRwdHTF+/HgUFRUZcS2IiIjIlAladoqLixEUFIR169bVOH358uVYs2YNNm7ciISEBNjb26N///4oLS3VzhMeHo5Lly4hNjYWhw4dwokTJzBx4kRjrQIRERGZOIlGo9EIHQIAJBIJ9u/fjyFDhgB4uFdHoVBgxowZmDlzJgBAqVTC1dUV27dvx8iRI5GWlob27dsjMTERwcHBAIDDhw/jlVdewZ07d6BQKOr03SqVCnK5HEqlEjKZzCDrR0Q1KCsCLn4F5N0A3DsBAX8BpFZCpyIiM1HX7belETPpJSMjA9nZ2QgLC9OOyeVy9OzZE/Hx8Rg5ciTi4+Ph6OioLToAEBYWBgsLCyQkJOC1116rcdllZWUoKyvTvlepVIZbESKq2f3rwI5BgOru72NuHYGIbwG7ZsLlIiLRMdkTlLOzswEArq6uOuOurq7aadnZ2XBxcdGZbmlpCScnJ+08NVmyZAnkcrn25eHh0cDpieiJjvxLt+gAQHYqcDJGmDxEJFomW3YMae7cuVAqldpXZmam0JGIGpeqCiD9+5qn/fIf42YhItEz2bLj5uYGAMjJydEZz8nJ0U5zc3NDbm6uzvTKykrk5eVp56mJjY0NZDKZzouIjEhiAVjUcm6O1Ma4WYhI9Ey27Pj4+MDNzQ1xcXHaMZVKhYSEBISEhAAAQkJCUFBQgKSkJO08R48ehVqtRs+ePY2emYjqyEIKdHi95mmdhhk3CxGJnqAnKBcVFeHatWva9xkZGUhJSYGTkxM8PT0RFRWFDz74AP7+/vDx8UF0dDQUCoX2iq2AgAAMGDAAEyZMwMaNG1FRUYHJkydj5MiRdb4Si4gE0n8x8Fs6cPfc72MBg4CQfwqXiYhEqd6Xnl+7dg3Xr1/Hc889Bzs7O2g0GkgkEr2W8eOPP+KFF16oNh4REYHt27dDo9Fg/vz52LRpEwoKCtCnTx+sX78ebdq00c6bl5eHyZMn49tvv4WFhQWGDh2KNWvWoGnTpnXOwUvPiQR0K/73S8/dOgqdhojMSF2333qXnfv372PEiBE4evQoJBIJ0tPT4evri8jISDRr1gwxMeZ3JQXLDhERkfmp6/Zb73N2pk2bBktLS9y+fRtNmjTRjo8YMQKHDx+uX1oiIiIiA9H7nJ3vv/8eR44cQatWrXTG/f39cevWrQYLRkRkNopygaTtQO5loEU7oNtYwKH2K0KJyLj0LjvFxcU6e3QeycvLg40NLxklMpqyQkB1D5B7ANbV/5skI/ntGrBtAFD86+9jZzcB474DWrQVLhcRael9GOvZZ5/Fzp07te8lEgnUajWWL19e48nGRNTA1Grg+2hgRRtgXQ8gpi3w4zKhUzVexz7QLToAUHIfiFsgTB4iqkbvPTvLly9HaGgozp07h/LycsyePRuXLl1CXl4eTp8+bYiMRPRHp1cDP635/X2ZCvhxMdC0BRAcKVisRuv6Uf3Gicjo9N6z06FDB1y9ehV9+vTB4MGDUVxcjNdffx3Jycnw8/MzREYi+qPELTWPn61lnAzL1rHmcT7MlMhk1OumgnK5HO+++25DZyGiuijKqWW89offkgF1HQMcXVjzOBGZBL3LzokTJx47/bnnnqt3GCKqA88Q4ObJ6uNevYyfhYA+0wDlHSD5M0BdCUikQOe/Ac/OFDoZEf2P3jcVtLCofuTrj3dOrqqqevpURsabCpJZyUwEdv4FqCj5fcxGDkR+B7gGCpersSvMfvj4C+fWgMxd6DREjUJdt99679nJz8/XeV9RUYHk5GRER0dj0aJF+iclIv14dAcm/ggkbHy4cXVpDzzzFuDkK3Syxs3BjffWITJRepcduVxebezFF1+EtbU1pk+frvMEciIykBZtgYGrhE5BRGQW9L4aqzaurq64cuVKQy2OiIiIqEHovWfn559/1nmv0WiQlZWFpUuXonPnzg2Vi4iIiKhB6F12OnfuDIlEgj+f1/zMM8/g008/bbBgRERERA1B77KTkZGh897CwgItWrSAra1tg4UiIiIiaih6lx0vLy9D5CAiIiIyiDqVnTVr1jx5pv+ZMmVKvcMQERERNbQ63VTQx8enbguTSHDjxo2nDmVsvKkgERGR+WnQmwr++TwdIiIiInPRYPfZISIiIjJF9Xrq+Z07d3Dw4EHcvn0b5eXlOtNWrlzZIMGIiIiIGoLeZScuLg5/+ctf4Ovri19++QUdOnTAzZs3odFo0LVrV0NkJCIiIqo3vQ9jzZ07FzNnzkRqaipsbW3x1VdfITMzE88//zyGDRtmiIxERERE9aZ32UlLS8OYMWMAAJaWlnjw4AGaNm2KBQsWYNmyZQ0ekIiIiOhp6F127O3ttefpuLu74/r169ppv/32W8MlIyIiImoAep+z88wzz+DUqVMICAjAK6+8ghkzZiA1NRVff/01nnnmGUNkJCIiIqo3vcvOypUrUVRUBAB4//33UVRUhC+++AL+/v68EouIiIhMjt5lZ/HixXjjjTcAPDyktXHjxgYPRURERNRQ9D5n59dff8WAAQPg4eGBWbNm4cKFC4bIRURERNQg9C4733zzDbKyshAdHY3ExER07doVgYGBWLx4MW7evGmAiERERI1IwW3gZAwQtwC4nSB0GlGo04NAH+fOnTv4/PPP8emnnyI9PR2VlZUNlc1o+CDQeqp4AFjZCZ2CiEg8Ln4NfD0RUFf8PhYcCQxcJVwmE1bX7fdTPRuroqIC586dQ0JCAm7evAlXV9enWRyZi3OfAqs6AovcgI+CgORdQiciIjJ/5cXAt1G6RQd4+DP3xnFBIolFvcrOsWPHMGHCBLi6umLs2LGQyWQ4dOgQ7ty509D5yNQk7wIOTQOUtx++z78JfDMJSP1S0FhERGbv5imgTFnztF/+Y9wsIqP31VgtW7ZEXl4eBgwYgE2bNmHQoEGwsbExRDYyRT99XMv4GqDjX42bhYhITKTWtU+zfMw0U1JVCZxeDZzfCZQWAH79gH7RgLOfoLH0Ljvvvfcehg0bBkdHRwPEIZOXl1HL+E2jxiAiEh3vPoCDO1CYVX1aRzN59uR/ZwBJ239/f2n/wz1Wb/8ENHURLJbeh7EmTJjAotOYuQfVPK6oZZyIiOpGagUM3wnYt/jDmDXQf0ntP3tNiSoLOP9Z9fHiX4GkHcbP8wd679mhRq7vHGD3cEBT9fuYhRXw3GzhMhERiYVHD2DaJSA99uEJy34vCLpHRC/303W3DX/0a5pxs/wJyw7pp3UYMOYb4PRHwG9XAJf2QJ9pgCefi0ZE1CAsbYCAgUKn0J9za0AirbnwNG9r/Dx/wLJD+vN59uGLiIjoEZkC6Dyq+u1ImjgD3cYKEukRlh0iIiJqGANXA7JWQPJnwIP8h1djhc4DHIS9D99T30FZDHgHZSIiIvNjlDsoExEREZk6lh0iIiISNZYdIiIiEjWWHSIiIhI1lh0iIiISNZYdIiIiEjWWHSIiIhI1lh0iIiISNZMuO1VVVYiOjoaPjw/s7Ozg5+eHhQsX4o/3QdRoNJg3bx7c3d1hZ2eHsLAwpKenC5iaiIgEoVY/fLr29oHApwOAn9YCFaVCpyITYNKPi1i2bBk2bNiAHTt2IDAwEOfOncO4ceMgl8sxZcoUAMDy5cuxZs0a7NixAz4+PoiOjkb//v1x+fJl2NraCrwGRERkNAcnAym7f39/Ox5I//7hw4slEuFykeBM+nERAwcOhKurK7Zu3aodGzp0KOzs7LBr1y5oNBooFArMmDEDM2fOBAAolUq4urpi+/btGDlyZI3LLSsrQ1lZmfa9SqWCh4cHHxdBRGSucn8B1veseVr4l4D/i8bNQ0YhisdF9OrVC3Fxcbh69SoA4MKFCzh16hRefvllAEBGRgays7MRFham/YxcLkfPnj0RHx9f63KXLFkCuVyufXl4eBh2RYiIyLAyE+o3jRoFkz6M9c4770ClUqFdu3aQSqWoqqrCokWLEB4eDgDIzs4GALi66j5N1dXVVTutJnPnzsX06dO17x/t2SEiIjMla/mYaQrj5SCTZNJlZ9++fdi9ezf27NmDwMBApKSkICoqCgqFAhEREfVero2NDWxsbBowqR7uXwe+jwbSjwBWTYCgkUDofMCmqTB5iIjEwK8f4OwP3P/TBSpNmgMd/ipMJjIZJn0Ya9asWXjnnXcwcuRIdOzYEaNHj8a0adOwZMkSAICbmxsAICcnR+dzOTk52mkm5UEBsO0V4Mp/AHUlUKYCzm4C/l3/4kZERAAsLIDRXwO+LwD438nILYOB0fsBW56L2diZ9J6dkpISWFjo9jGpVAq1Wg0A8PHxgZubG+Li4tC5c2cADw9JJSQk4O233zZ23Cf7+QugqIbDa9d+ALJTAbeOxs9ERCQWjp7AmANA8X1AXQE4mOAvvSQIky47gwYNwqJFi+Dp6YnAwEAkJydj5cqViIyMBABIJBJERUXhgw8+gL+/v/bSc4VCgSFDhggbvib3rz1+GssOEdHTs3cWOgGZGJMuOx9//DGio6Pxj3/8A7m5uVAoFPj73/+OefPmaeeZPXs2iouLMXHiRBQUFKBPnz44fPiwad5jxzXwMdM6GC8HERFRI2LS99kxlrpep//UyouBjX2AvBu64x2GAn/91HDfS0RUm/ybQPZFwMkXcG0vdBoivdR1+23Se3ZEx9oeGPcd8ONS4Orh/12NNQroEyV0MiJqbNRVwMEp/7vj8P9+520dBgzbDtg4CJmMqMGx7BibgxswaLXQKYiosTuzHkjZpTt27Qcgdh4wcJUwmYgMxKQvPSciIgNJ+bzm8QtfADy7gUSGZYeIqDEqL6x5vKLk4SEuIhFh2SEiaoz8+9c83joUkPIMBxIXlh0iosbo+dmAc2vdsSbNgZc+ECYPkQGxvhMRNUZNXYC/nwB+3gdkXXh46XnncN6Qj0SJZYeIqLGytgeCxwmdgsjgeBiLiIiIRI1lh4iIiESNZYeIiIhEjWWHiIiIRI0nKJP4lRUBad8CD/IAn+cBNz5hnoioMWHZIXG7kwTs/uvDovNIcCSf/UNE1IjwMBaJl0YDfD1Bt+gAwLlPgbRDwmQiIiKjY9kh8cpOBfKu1zzt0n7jZiEiIsGw7JCIPe7JzXyqMxFRY8GyQ+Ll1unhLfBrEviacbMQEZFgWHZIvCQS4LVNgF0z3fGuEUC7gcJkIiIio+PVWCRuHt2BqFTg8kGg5D7g2xdw7yR0KiIiMiKWHRI/GwegS7jQKYiISCA8jEVERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESiZvJl5+7du3jjjTfg7OwMOzs7dOzYEefOndNO12g0mDdvHtzd3WFnZ4ewsDCkp6cLmJiI6uJuwQPEfH8FU/cmY/OJG1A+qBA6EhGJlEmXnfz8fPTu3RtWVlb47rvvcPnyZcTExKBZs2baeZYvX441a9Zg48aNSEhIgL29Pfr374/S0lIBkxPR46RkFqD/qhP4+Og1fJNyD4v+m4aBH59Eror/3RJRw5NoNBqN0CFq88477+D06dM4efJkjdM1Gg0UCgVmzJiBmTNnAgCUSiVcXV2xfft2jBw5ssbPlZWVoaysTPtepVLBw8MDSqUSMpms4VeEiHQM2/gTEm/mVxuPCPHC+4M7CJCIiMyRSqWCXC5/4vbbpPfsHDx4EMHBwRg2bBhcXFzQpUsXbN68WTs9IyMD2dnZCAsL047J5XL07NkT8fHxtS53yZIlkMvl2peHh4dB14OIfvegvKrGogMAJ9J/M3IaImoMTLrs3LhxAxs2bIC/vz+OHDmCt99+G1OmTMGOHTsAANnZ2QAAV1dXnc+5urpqp9Vk7ty5UCqV2ldmZqbhVoKIdFhJJWhiLa1xmszW0shpiKgxMOmfLGq1GsHBwVi8eDEAoEuXLrh48SI2btyIiIiIei/XxsYGNjY2DRWTiPRgKbXA611bYteZ29WmDe/OvaxE1PBMes+Ou7s72rdvrzMWEBCA27cf/pB0c3MDAOTk5OjMk5OTo51GRKbnX68EYECgGySSh++tpRaY8KwP/tbDU9hgRCRKJr1np3fv3rhy5YrO2NWrV+Hl5QUA8PHxgZubG+Li4tC5c2cAD09WSkhIwNtvv23suERUR02sLbFxdDfcvl+CzPwStHVzQPOm3NtKRIZh0mVn2rRp6NWrFxYvXozhw4fj7Nmz2LRpEzZt2gQAkEgkiIqKwgcffAB/f3/4+PggOjoaCoUCQ4YMETY8ET2Rp3MTeDo3EToGEYmcSZed7t27Y//+/Zg7dy4WLFgAHx8frF69GuHh4dp5Zs+ejeLiYkycOBEFBQXo06cPDh8+DFtbWwGTExERkakw6fvsGEtdr9MnIiIi0yGK++wQERERPS2WHSIiIhI1lh0iIiISNZYdIiIiEjWWHSIiIhI1k770nIjI1JxK/w07428iR1WKYG8nTHjWF25y3uqCyJSx7BAR1dG+c5mY/eXP2vcX7ijxn5+z8M3k3nCVsfAQmSoexiIiqoPKKjVWHLlSbTxbVYpPT2UIkIiI6oplh4ioDu7kP0BuYVmN087fzjdyGiLSB8sOEVEdODW1hrW05h+ZbnI7I6chIn2w7BAR1YHM1gqvd21ZbdxCAkSEeAmQiIjqiicoExHV0Xt/CYSFhQRfJd1BWaUarZrZYc6Adgj2dhI6GhE9Bh8ECj4IlIj0U1xWCeWDCrjJbGFhIRE6DlGjVdftN/fsEBHpyd7GEvY2/PFJZC54zg4RERGJGssOERERiRrLDhEREYkayw4RERGJGssOkZmpqFKDF1ESEdUdLycgMhPJt/Ox5L+/4OzNPDjYWmJkdw/M7N8WNpZSoaMREZk0lh0iM3D7fgne2JKA4vIqAEBhaSU2n8xAXnEFYoYHCZyOiMi08TAWkRnYnXBLW3T+6EDKXeSqSgVIRERkPlh2iMzAzfvFNY5XqTXIzC8xchoiIvPCskNkBgIV8hrHba0s0LqFg5HTEBGZF5YdIjPwt56ecJXZVBsf19sH8iZWAiQiIjIfPEGZyAw0b2qDL9/qhbVHr+HUtd/gZG+NkT08EN7TS+hoREQmj2WHyEx4ODXBsr92EjoGEZHZ4WEsIiIiEjXu2SEiIhJQRZUaXyRm4silbFhaSDC4c0sM7qyARCIROpposOwQEREJRKPR4O+fJeHoL7nasWNXfkVCxn0seZ2HrRsKD2MREREJ5PS1+zpF55HPz2biWm6RAInEiWWHiIhIIGdv5tU6LfEx00g/LDtEREQCcXGofv+sukwj/bDsEBERCWRQkAKONdwY1NOpCZ5v00KAROLEE5SpXkrKK/FZ/C0c/SUXTayleL1rKwwKUggdi4jIrMjtrLAzsgfmfJWKtCwVAKC7dzMs/2sQLKXmtz+itKIK313MQrayDMHezdDd20noSABYdqgeyivVeGNLAs7fLtCOHbvyKy7eU2LuywHCBSMiMkOdWjniu6nP4vb9ElhKJVA42gkdqV7ScwoRviUBuYVl2rEX27tifXhXWAlc3MyvNpLgvruYpVN0Htl6MgM5qlLjByIiEgFP5yZmW3QAYM5XP+sUHQCIvZyDvWdvC5Todyw7pLfarhCoVGuQXEMJIiIicctRldb4SzAA/Cc1y7hhasCyQ3pzk9nWOq2mJ3MTEZG4Pe5ezxYmcCdolh3S21+7ecDeWlptvFMrObp4NhMgERERCclFZovu3jX//H+1k7uR01THskN6c5PbYtu4Hmjr6gAAkEiAF9q2wJYxwQInIyIioSwd2gkt/3TO0aud3DGyu6dAiX4n0Wg0GqFDCE2lUkEul0OpVEImkwkdx6xk5pXAzlqK5k15+IqIqLErr1Qj9nIOslWlCPZqhiAPR4N+X12337z0nJ6Kh1MToSMQEZGJsLa0MInDVn/Gw1hERCQ6uapSXLqnRFllldBRyARwzw4REYlGYWkFZn/5M45cyoZaAzjZW2PmS23xt57CnzdCwuGeHSIiEo25X6fiu4sPiw4A5BWX41/7U/HTtd+EDUaCYtkxgvzicqTeUUJVWiF0FCIi0corLsd3F7NrnLbbBO7iS8LhYSwDqqxSY8Ghy9h7NhPlVWrYWllgbC8fzBnQFhITuMkSEZGY5JeUo0pd8wXG94vKahynxsGs9uwsXboUEokEUVFR2rHS0lJMmjQJzs7OaNq0KYYOHYqcnBzhQv7B2mPXsDP+Fsqr1ACA0go1Nh6/jp3xtwRORkQkPt7O9nCX13yH92d8nY2chkyJ2ZSdxMREfPLJJ+jUqZPO+LRp0/Dtt9/i3//+N44fP4579+7h9ddfFyilrt0JNe823Z3AskNE1NCkFhK8+2oApBa6e859m9tjbC9vYUKRSTCLw1hFRUUIDw/H5s2b8cEHH2jHlUoltm7dij179qBfv34AgG3btiEgIABnzpzBM888I1RkAA+PH9fkflHN40RE9HQGdlLAo1kT7Em4jZzCUnT3dsIbPb0gb2IldDQSkFmUnUmTJuHVV19FWFiYTtlJSkpCRUUFwsLCtGPt2rWDp6cn4uPjay07ZWVlKCv7/fitSqUySO4QX2ecquEKgBA/7k4lIjKUIA9Hg9+5l8yLyR/G2rt3L86fP48lS5ZUm5adnQ1ra2s4OjrqjLu6uiI7u+Yz8gFgyZIlkMvl2peHh0dDxwYAzB7QFk1tdPtksyZWiAprY5DvIyIioupMes9OZmYmpk6ditjYWNja1nzSWX3MnTsX06dP175XqVQGKTydWjniP1P6YGf8Ldz4tQjt3GUYE+IFd7ndkz9MRGQAWcoH+CIxE3fzH6CLZzO81qUl7KylQsciMiiTLjtJSUnIzc1F165dtWNVVVU4ceIE1q5diyNHjqC8vBwFBQU6e3dycnLg5uZW63JtbGxgY2OcB1d6OdsjemB7o3wXEdHjJN3KQ8SniSgqqwQA/DvpDrb/lIF9fw+BYxNrgdMRGY5JH8YKDQ1FamoqUlJStK/g4GCEh4dr/2xlZYW4uDjtZ65cuYLbt28jJCREwORERKbnvYOXtUXnkas5RdhyMkOgRETGYdJ7dhwcHNChQwedMXt7ezg7O2vHx48fj+nTp8PJyQkymQz//Oc/ERISIviVWEREpuR+URlS7yprnHbsSi5m9m9r5ERExmPSZacuVq1aBQsLCwwdOhRlZWXo378/1q9fL3QsIiKTYmslhZVUgoqq6ncYdrA1+00B0WNJNBpNzffWbkRUKhXkcjmUSiVkMpnQcaiB5ReXY3/yXfxaVIYe3k54vk0LWFjwcR3U+Ez5PBkHL9yrNh4zLAhDu7USIBHR06nr9pt1nkTt/O18RHx6FoWlD89T2IDreK5NC2wZEwxrS5M+ZY2owS0c3AF5xeXa+39ZSSWI7O3DokOix7JDovbOVz9ri84jJ67+ii/OZWL0M14CpSIShryJFXa92RNXcwpxt+ABAhUyuDg03G09iEwVf7Ul0br5WzGu5hTVOO37S7XfdJJI7Nq4OuCFti4sOtRosOyQaD3uMJUND2ERETUa/IlPoqVwtEMPb6capw3u3NLIaYiISCgsOyRqK4YFwbe5vfa9RAKM7eWNQUEKAVMREZEx8QRlEjVP5yb4YfrzOHXtN+QWlqG7dzN4Ods/+YNERCQaLDskehYWEjzXpoXQMYiISCA8jEVERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosayQ0RERKLGskNERESixrJDREREosZnYwHQaDQAAJVKJXASIiIiqqtH2+1H2/HasOwAKCwsBAB4eHgInISIiIj0VVhYCLlcXut0ieZJdagRUKvVuHfvHhwcHCCRSBp8+SqVCh4eHsjMzIRMJmvw5RsL18N0iGEdAHGshxjWAeB6mBIxrANgnPXQaDQoLCyEQqGAhUXtZ+Zwzw4ACwsLtGrVyuDfI5PJzPr/uI9wPUyHGNYBEMd6iGEdAK6HKRHDOgCGX4/H7dF5hCcoExERkaix7BAREZGosewYgY2NDebPnw8bGxuhozwVrofpEMM6AOJYDzGsA8D1MCViWAfAtNaDJygTERGRqHHPDhEREYkayw4RERGJGssOERERiRrLDhEREYkay46BnThxAoMGDYJCoYBEIsGBAweEjqS3JUuWoHv37nBwcICLiwuGDBmCK1euCB1LLxs2bECnTp20N7cKCQnBd999J3Ssp7Z06VJIJBJERUUJHUUv7733HiQSic6rXbt2QsfS2927d/HGG2/A2dkZdnZ26NixI86dOyd0LL14e3tX+7eQSCSYNGmS0NHqrKqqCtHR0fDx8YGdnR38/PywcOHCJz4vyRQVFhYiKioKXl5esLOzQ69evZCYmCh0rMd60nZOo9Fg3rx5cHd3h52dHcLCwpCenm7UjCw7BlZcXIygoCCsW7dO6Cj1dvz4cUyaNAlnzpxBbGwsKioq8NJLL6G4uFjoaHXWqlUrLF26FElJSTh37hz69euHwYMH49KlS0JHq7fExER88skn6NSpk9BR6iUwMBBZWVna16lTp4SOpJf8/Hz07t0bVlZW+O6773D58mXExMSgWbNmQkfTS2Jios6/Q2xsLABg2LBhAieru2XLlmHDhg1Yu3Yt0tLSsGzZMixfvhwff/yx0NH09uabbyI2NhafffYZUlNT8dJLLyEsLAx3794VOlqtnrSdW758OdasWYONGzciISEB9vb26N+/P0pLS40XUkNGA0Czf/9+oWM8tdzcXA0AzfHjx4WO8lSaNWum2bJli9Ax6qWwsFDj7++viY2N1Tz//POaqVOnCh1JL/Pnz9cEBQUJHeOpzJkzR9OnTx+hYzS4qVOnavz8/DRqtVroKHX26quvaiIjI3XGXn/9dU14eLhAieqnpKREI5VKNYcOHdIZ79q1q+bdd98VKJV+/rydU6vVGjc3N82HH36oHSsoKNDY2NhoPv/8c6Pl4p4d0ptSqQQAODk5CZykfqqqqrB3714UFxcjJCRE6Dj1MmnSJLz66qsICwsTOkq9paenQ6FQwNfXF+Hh4bh9+7bQkfRy8OBBBAcHY9iwYXBxcUGXLl2wefNmoWM9lfLycuzatQuRkZEGeSiyofTq1QtxcXG4evUqAODChQs4deoUXn75ZYGT6aeyshJVVVWwtbXVGbezszO7PZ+PZGRkIDs7W+dnlVwuR8+ePREfH2+0HHwQKOlFrVYjKioKvXv3RocOHYSOo5fU1FSEhISgtLQUTZs2xf79+9G+fXuhY+lt7969OH/+vMkfx3+cnj17Yvv27Wjbti2ysrLw/vvv49lnn8XFixfh4OAgdLw6uXHjBjZs2IDp06fjX//6FxITEzFlyhRYW1sjIiJC6Hj1cuDAARQUFGDs2LFCR9HLO++8A5VKhXbt2kEqlaKqqgqLFi1CeHi40NH04uDggJCQECxcuBABAQFwdXXF559/jvj4eLRu3VroePWSnZ0NAHB1ddUZd3V11U4zBpYd0sukSZNw8eJFs/wto23btkhJSYFSqcSXX36JiIgIHD9+3KwKT2ZmJqZOnYrY2Nhqv/2Zkz/+xt2pUyf07NkTXl5e2LdvH8aPHy9gsrpTq9UIDg7G4sWLAQBdunTBxYsXsXHjRrMtO1u3bsXLL78MhUIhdBS97Nu3D7t378aePXsQGBiIlJQUREVFQaFQmN2/xWeffYbIyEi0bNkSUqkUXbt2xahRo5CUlCR0NLPGw1hUZ5MnT8ahQ4dw7NgxtGrVSug4erO2tkbr1q3RrVs3LFmyBEFBQfjoo4+EjqWXpKQk5ObmomvXrrC0tISlpSWOHz+ONWvWwNLSElVVVUJHrBdHR0e0adMG165dEzpKnbm7u1crygEBAWZ3OO6RW7du4YcffsCbb74pdBS9zZo1C++88w5GjhyJjh07YvTo0Zg2bRqWLFkidDS9+fn54fjx4ygqKkJmZibOnj2LiooK+Pr6Ch2tXtzc3AAAOTk5OuM5OTnaacbAskNPpNFoMHnyZOzfvx9Hjx6Fj4+P0JEahFqtRllZmdAx9BIaGorU1FSkpKRoX8HBwQgPD0dKSgqkUqnQEeulqKgI169fh7u7u9BR6qx3797VbsFw9epVeHl5CZTo6Wzbtg0uLi549dVXhY6it5KSElhY6G7OpFIp1Gq1QImenr29Pdzd3ZGfn48jR45g8ODBQkeqFx8fH7i5uSEuLk47plKpkJCQYNRzJnkYy8CKiop0flvNyMhASkoKnJyc4OnpKWCyups0aRL27NmDb775Bg4ODtrjrHK5HHZ2dgKnq5u5c+fi5ZdfhqenJwoLC7Fnzx78+OOPOHLkiNDR9OLg4FDtXCl7e3s4Ozub1TlUM2fOxKBBg+Dl5YV79+5h/vz5kEqlGDVqlNDR6mzatGno1asXFi9ejOHDh+Ps2bPYtGkTNm3aJHQ0vanVamzbtg0RERGwtDS/zcKgQYOwaNEieHp6IjAwEMnJyVi5ciUiIyOFjqa3I0eOQKPRoG3btrh27RpmzZqFdu3aYdy4cUJHq9WTtnNRUVH44IMP4O/vDx8fH0RHR0OhUGDIkCHGC2m0674aqWPHjmkAVHtFREQIHa3OasoPQLNt2zaho9VZZGSkxsvLS2Ntba1p0aKFJjQ0VPP9998LHatBmOOl5yNGjNC4u7trrK2tNS1bttSMGDFCc+3aNaFj6e3bb7/VdOjQQWNjY6Np166dZtOmTUJHqpcjR45oAGiuXLkidJR6UalUmqlTp2o8PT01tra2Gl9fX827776rKSsrEzqa3r744guNr6+vxtraWuPm5qaZNGmSpqCgQOhYj/Wk7ZxardZER0drXF1dNTY2NprQ0FCj/39NotGY4S0miYiIiOqI5+wQERGRqLHsEBERkaix7BAREZGosewQERGRqLHsEBERkaix7BAREZGosewQERGRqLHsEBERkaix7BCRyenbty+ioqKEjkFEIsGyQ0RERKLGskNERESixrJDRCZJrVZj9uzZcHJygpubG9577z3ttNu3b2Pw4MFo2rQpZDIZhg8fjpycHO30sWPHVnuiclRUFPr27at9/+WXX6Jjx46ws7ODs7MzwsLCUFxcrJ2+ZcsWBAQEwNbWFu3atcP69esNtapEZGAsO0Rkknbs2AF7e3skJCRg+fLlWLBgAWJjY6FWqzF48GDk5eXh+PHjiI2NxY0bNzBixIg6LzsrKwujRo1CZGQk0tLS8OOPP+L111/Ho+ci7969G/PmzcOiRYuQlpaGxYsXIzo6Gjt27DDU6hKRAVkKHYCIqCadOnXC/PnzAQD+/v5Yu3Yt4uLiAACpqanIyMiAh4cHAGDnzp0IDAxEYmIiunfv/sRlZ2VlobKyEq+//jq8vLwAAB07dtROnz9/PmJiYvD6668DAHx8fHD58mV88skniIiIaND1JCLD454dIjJJnTp10nnv7u6O3NxcpKWlwcPDQ1t0AKB9+/ZwdHREWlpanZYdFBSE0NBQdOzYEcOGDcPmzZuRn58PACguLsb169cxfvx4NG3aVPv64IMPcP369YZbQSIyGu7ZISKTZGVlpfNeIpFArVbX6bMWFhbaQ1KPVFRUaP8slUoRGxuLn376Cd9//z0+/vhjvPvuu0hISECTJk0AAJs3b0bPnj11liGVSuuzKkQkMO7ZISKzEhAQgMzMTGRmZmrHLl++jIKCArRv3x4A0KJFC2RlZel8LiUlRee9RCJB79698f777yM5ORnW1tbYv38/XF1doVAocOPGDbRu3Vrn5ePjY/D1I6KGxz07RGRWwsLC0LFjR4SHh2P16tWorKzEP/7xDzz//PMIDg4GAPTr1w8ffvghdu7ciZCQEOzatQsXL15Ely5dAAAJCQmIi4vDSy+9BBcXFyQkJODXX39FQEAAAOD999/HlClTIJfLMWDAAJSVleHcuXPIz8/H9OnTBVt3Iqof7tkhIrMikUjwzTffoFmzZnjuuecQFhYGX19ffPHFF9p5+vfvj+joaMyePRvdu3dHYWEhxowZo50uk8lw4sQJvPLKK2jTpg3+7//+DzExMXj55ZcBAG+++Sa2bNmCbdu2oWPHjnj++eexfft27tkhMlMSzZ8PbBMRERGJCPfsEBERkaix7BAREZGosewQERGRqLHsEBERkaix7BAREZGosewQERGRqLHsEBERkaix7BAREZGosewQERGRqLHsEBERkaix7BAREZGo/T9WL38VBPwWowAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "strip = sns.stripplot(data=df, x=\"house\", y=\"value\", hue=\"group\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "结果变量为大鼠 MT 含量（数值变量），影响变量有区组和不同的处理（2 $ \\times $ 多项无序分类变量）。统计时应剔除区组的影响，分析不同处理对 MT 的含量的影响。\n",
    "\n",
    "可以使用 Two way ANOVA, 分别推测两种影响变量对结果变量的影响。\n",
    "\n",
    "### 假设\n",
    "\n",
    "$ H_0 $: $ \\mu _1 = \\mu _2 = \\mu _3 $，即三组大鼠的 MT 含量的总体均值相同。  \n",
    "$ H_1 $: $ \\exists i,j \\in \\{1,2,3\\}: \\mu _i \\neq \\mu _j $ , 即三组大鼠的 MT 含量的总体均值不全相同。\n",
    "\n",
    "### 检验"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sum_sq</th>\n",
       "      <th>df</th>\n",
       "      <th>F</th>\n",
       "      <th>PR(&gt;F)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>C(group)</th>\n",
       "      <td>33078.798</td>\n",
       "      <td>2.0</td>\n",
       "      <td>341.918569</td>\n",
       "      <td>4.800931e-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C(house)</th>\n",
       "      <td>1276.963</td>\n",
       "      <td>9.0</td>\n",
       "      <td>2.933180</td>\n",
       "      <td>2.485531e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Residual</th>\n",
       "      <td>870.702</td>\n",
       "      <td>18.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             sum_sq    df           F        PR(>F)\n",
       "C(group)  33078.798   2.0  341.918569  4.800931e-15\n",
       "C(house)   1276.963   9.0    2.933180  2.485531e-02\n",
       "Residual    870.702  18.0         NaN           NaN"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "formula = \"value ~ C(group) + C(house)\"\n",
    "\n",
    "model = ols(formula, data=df).fit()\n",
    "\n",
    "table = sm.stats.anova_lm(model, typ=2)\n",
    "\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于 C(group)，PR(>F) —— 即 P-value —— 小于 0.05，差别具有统计学意义，因而拒绝 $ H_0 $，认为三组大鼠的 MT 含量的总体均值不全相同，即不同处理对大鼠 MT 含量有影响。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
